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1 Introduction 



The magnetic induction equations are a special form of the Maxwell's equations that 
describe the evolution of the magnetic field under the influence of a given velocity 
field. These equations arise in a wide variety of applications in plasma physics, as- 
trophysics and electrical engineering. One important application are the equations of 
magneto-hydro dynamics (MHD). These equations combine the Euler equations of 
gas dynamics with the magnetic induction equations. Our goal in this paper is to de- 
scribe stable and high-order accurate numerical schemes for the magnetic induction 
equations. 

We start with a brief description of how the equations are derived. Let the mag- 
netic field and given velocity field be denoted by B and u respectively. Faraday's law 
for the magnetic flux across a surface S bounded by a curve dSis given by (see [ ]), 



^ J B-dS = ^E-dl. 



dt . 

s ds 

Using the Stokes theorem and the fact that the electric field, E, in a co-moving frame 
is zero and the magnetic resistivity is zero, Faraday's law takes the form, 

^+div(u(g)B-B(g)u) = -udiv(B). (1.1) 
ot 

Using simple vector identities, (1.1) can be rewritten as, 

5,B + curl(B X u) = -udiv(B). (1.2) 



Magnetic monopoles have never been observed in nature. As a consequence, the mag- 
netic field is always assumed to be divergence free, i.e., div(B) = 0. Hence, it is com- 
mon to set the right-hand side of (1.2) to zero and couple the induction equation with 
the divergence constraint in order to obtain 

5,B + curl(Bxu) =0, 

^ ^ (1.3) 
div(B) =0,B(x,0) =Bo(x). 

This form (1.3) is commonly used in the literature as the appropriate form of the 
magnetic induction equations to study and discretize. It is easy to see that (1.3) is 
hyperbolic but not strictly hyperbolic. An important tool in the analysis of hyper- 
bolic system of equations is the derivation of energy estimates. The usual procedure 
in deriving energy estimates consists of symmetrizing the hyperbolic system. It is 
not possible to symmetrize (1.3) without explicitly using the divergence constraint. 
Hence, it is difficult to obtain energy stability starting from (1.3). 
On the other hand, we can use the following vector identity 



curl(B X u) = Bdivu - udiv(B) + (u • V) B - (B • V) u 

= {u'B)^ + {u^B)^ + (^^B)^ - udiv(B) - (B • V)u, 
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and rewrite (1.1) in the non-conservative symmetric form, 



a,B + (u-V)B 



= -B(divu) + (B-V)u 
= M(Z)u)B, 



(1.4) 



where the Du denotes the gradient of u and the matrix M{Du) is given by 



( 



dyU^ — d^u^ 



dyU^ 

dyU^ 



dxU^ — dyU^ 



M{Du) = 




Introducing the matrix, 



C= - 



djU^ dyU^ d^u^ 
djcU^ dyU^ d^u^ 
dxU^ dyU^ dzU^ 



) 



(1.1) can also be written in the following "conservative" symmetric form. 



dt^ + dx (A^B) + dy (A^B) + (A^B) + CB = 0, 



(1.5) 



where A^ = w7 for / = 1 , 2, 3. Note that the symmetrized matrices in (1 .5) are diagonal 
and that the only coupling in the equations is through the lower order terms. These 
symmetrized forms are in the same spirit as the non-linear symmetrized forms of 
MHD equations introduced in [ ] . 

Furthermore, by taking divergence on both sides of (1.2) we get 



Hence, if div(Bo(x)) = 0, also div(B(x,^)) = for ^ > 0. This implies that all the 
above forms (1 .5), and (1.3) are equivalent (at least for smooth solutions). Introducing 
the space H^^^ as 



we have the following theorem: 

Theorem 1.1 Assume that the velocity field u is sufficiently smooth, and that Bq G 



7^^^^(R^). Then there exists a unique weak solution B G C([0, r];if^^^(R^)) of (1.5). 
The solution B satisfies the energy estimate, 



The constant Cj depends only on the final time T. Furthermore, //'div(Bo) = 0, then 
the physical form (1.1) and the symmetric form (1.5) are equivalent to the constrained 
form (1.3), i.e., B is also the unique weak solution of (1.3). 



(div(B)), +div(udiv(B)) = 0. 



(1.6) 



H^^{R^) = {w:R^ \ |w| GL^(M^), div ( w) G L^(M^)} , 



ll^(^^)ll//div(M3) < Ct ||Bo||^div(]^3) 
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The proof of the above theorem uses the energy estimate and we will provide a sketch 
of the proof for the two-dimensional version of the equations together with boundary 
conditions later in this paper. 

Even though the magnetic induction equations are linear, the presence of variable 
coefficients and lower order terms means that general closed form solutions are not 
available. Hence, one has to design suitable numerical schemes for these equations. 
Furthermore, since these equations appear as a sub-model in the MHD equations, 
the design of stable and high-order accurate numerical schemes for the induction 
equations can lead to the design of robust schemes for the non-linear MHD equations. 

Most of the attention in the literature has been focused on the constrained form 
(1 .3). The key issue in the design of a suitable numerical scheme to approximate (1 .3) 
has been the treatment of the divergence constraint. A widely used approach has been 
to employ projection methods based on a Hodge decomposition of the magnetic field. 
A base (finite difference or finite volume) scheme is used to evolve the magnetic 
field. The evolved field, which need not be divergence free, is then corrected for 
divergence errors by solving an elliptic equation (see [ ]). The resulting method is 
computationally expensive, as the elliptic equation has to be solved at every time 
step. 

Another common approach is to discretize (1.3) such that some particular form 
of discrete divergence is preserved at each time step (see [ ]). This approach is 
equivalent to staggering the velocity and magnetic fields in each direction (see [4, 
1,20,5] and a detailed comparison in [2 ]). Some of these schemes are proved to 
be von Neumann stable in the special case of constant velocity fields. No stability 
analysis is available either in the case of variable velocity fields or for problems with 
boundary conditions. These schemes also involve wider stencils than what is required 
for a standard finite difference scheme. 

Despite all the attempts at finding a suitable discretization of (1.3) and preserving 
a special form of discrete divergence, it is not clear as to whether such an approach is 
appropriate. Furthermore, there are many different choices for the discrete divergence 
operator and preserving some form of discrete divergence exactly does not lead to 
preservation or even keeping divergence errors small for a different form. The main 
aim should be to design a stable scheme to approximate magnetic fields and it is not 
clear whether preserving divergence in a particular discrete form helps. One reason 
for the difficulties in proving stability of discretizations for (1.3) with general velocity 
fields may lie in the very form of these equations. As remarked earlier, (1.3) are not 
symmetrizable directly and thus one cannot obtain energy estimates in this form. This 
remains true for discretizations of (1.3). 

A different approach consisting of discretizing the physical form (1.1) was pro- 
posed in [17] for the non-linear MHD equations. Adapting this to (1.1) implies using 
a standard upwind scheme for the convection part and a centered discretization of the 
source terms. From (1.6), one can expect that divergence errors will be transported 
out of the domain for transparent boundary conditions. This approach does not im- 
ply stability either and can lead to oscillations as reported in [ ]. A discontinuous 
Galerkin based discretization of the symmetric form (1.5) was proposed in [ ]. 

In a recent paper [ ], the authors discretized the symmetric form (1.4) by using 
a first order accurate upwind finite difference scheme. The resulting scheme also 
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implied an upwind discretization of the convection term in (1.1) with an upwind 
discretization of the source term. This scheme was shown to be energy stable even 
with variable velocity fields and to be TVD for constant velocity fields. 

Furthermore, boundary conditions were not considered either in [6] or any of 
the aforementioned papers. High-order accurate schemes will lead to much better 
resolution of interesting solution features and a stable discretization of the boundary 
conditions (while still preserving high order of accuracy) is desirable. 

Our aim in this paper is to design stable and high-order accurate schemes for 
initial-boundary value problems corresponding to the magnetic induction equations 
by discretizing the non-conservative symmetric form (1.4). The spatial derivatives are 
approximated by second and fourth-order SBP (Summation-By-Parts) operators. The 
boundary conditions are weakly imposed by using a SAT (Simultaneous Approxima- 
tion Term) and time integration is performed by standard Runge-Kutta schemes. The 
SBP-SAT framework has been used to obtain stable and accurate high order schemes 
for a wide variety of hyperbolic problems in recent years. See [ ] and the references 
therein for more details. 

The SBP-SAT schemes use centered finite difference stencils in the interior, which 
lead to oscillations in the vicinity of discontinuities. We apply well-known SBP-SAT 
compatible numerical diffusion operators in case of discontinuous data. 

The rest of this paper is organized as follows: In Section 2, we state the energy 
estimate for the initial-boundary value problem corresponding to (1.4) in order to 
motivate the proof of stability for the scheme. In Section 3, we present the SBP-SAT 
scheme and show stability. Numerical experiments are presented in Section 4 and 
conclusions are drawn in Section 5. 



2 The Continuous problem 

For ease of notation, we shall restrict ourselves to two spatial dimensions in the re- 
mainder of this paper. Extending the results to three dimensions is straightforward. 
In two dimensions, the non-conservative symmetric form (1.4) reads 



where 



B, +AiB;, +A2B. -CB = 0, (2.1) 



with B = {B^ ^B^Y and u = {u^ ^u^Y denoting the magnetic and velocity fields re- 
spectively. In component form, (2.1) becomes 

{B' \ + u' {B' + u\B' )y = -{u\B' + {u' )yB^ 
{B^)t + u' {b\ + u\b\ = {u\b' - {u\b\ 



To begin with, we shall consider (2.1) in the domain (x,};) G 12 = [0, 1]^. 
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We augment (2.1) with initial conditions, 



B(x,0)=Bo(x) xeQ, 



(2.3) 



and Dirichlet boundary conditions, 



(2.4) 



where 1a denotes the characteristic function of the set A. Note that we only impose 
boundary conditions on the set where the characteristics are entering the domain. 

Definition 2.1 Weak solution: A function B : 12 ^ such that B eC{[0,T];H^{n)) 
is defined as a weak solution of (2.1) with initial data (2.3) and boundary data (2.4) if 
it satisfies the weak formulation of (2.1) in Q, i.e.. 



for all test functions (jO G C^(I2 x [0, T)).By TrB we mean the trace of B at the 
boundary. The boundary conditions (2.4) are taken in the sense of traces. 

We shall always assume that the initial and boundary data satisfy the compatibility 
conditions, i.e., specific criteria that guarantee smoothness of the solution, see [ ]. 

Theorem 2.1 Assume that Bq G {Q), that geH^ {dQ x [0, T])forT >0 and that 
and are in H^{Q x [0, T]). Then there exists a function B G C([0, r],L^(I2)) fl 
U°{[OJ]\H^{Q.)) which is the unique weak solution of (2.1) with the initial and 
boundary conditions (2.3) and (2.4). 

Furthermore, it satisfies the following stability estimate 



where a is a positive constant. 

Proof The proof of this theorem is standard. Assume first that g, Bq and u are in 
C°°. Since the compatibility conditions are satisfied, a unique solution exists by the 




(2.5) 




(2.6) 
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method of characteristics. Let (aVO) = max{a,0} and {a AO) = inin{a,0}. Multi- 
plying the equation by B and integrating over Q yields 



d 
dt 



f B^Bdxdy 
Jn 

= J^B{2C^diy{u))Bdxdy- u^Tr{B^B) ^>^ + ^ i/^rr(B^B) dx 

<c [ B^Bdxdy 
Jn 

+ {u\0,y,t)VO) (rr(B^B)) dy- J^' {u\hy,t) AO) (rr(B^B)) dy 

^ {u^{x,0,t)VO) (rr(B^B)) dx- {u^{x,l,t) AO) (rr(B^B)) dx 

<c( [ {B^B)dxdy^ I g^ds 
\Jn Jdn ) 



for some constant c depending on u and its first derivatives. Via the Gronwall inequal- 
ity we get the bound 



2 



l|B(-,f)lll2(i2)<^'^'(^l|Bo||^2(^)+ j-^^rdsdt 

Set P = Bx and Q = B^, applying to (2.1) yields 

P, + j^^Px + u% = ulV + ulQ + CP + GB. (2.7) 

Furthermore, P{x^y^ 0) = dxBo{x^y) and at those parts of dQ where we impose bound- 
ary data 

u^F = Cg — gt — li^^y onx = Oandx=l, 

w^Q = Cg-gt- u^gx on J = and J = 1. 

We shall also be needing P on _y = and 1 and Q on x = and 1 . These are given by 
gx and g^; respectively. 

Multiplying (2.7) with 2P^ and rearranging yields 



+ (w^P^)^ + {u^^^)y = -^xP^ - 2i^xP^Q + 2P^CP + 2P^QB. 
We also have an analogous equation for Q^; 

Qj + (u^Q^)^ + (^^Q^)^ = -^yQ^ - 24p^Q + 2P^^Q + 2Q^C^B. 
Adding these two equations we find 

(P2 + Q2)^ + («1 (P2 + Q2) )^ + („2 (p2 + q2) )^ ^ (2.8) 

where by Holder's inequality R has the bound 



J^R{x,y,t)dxdy < c (||B(-,0||'2(^) + ||P(-,0IIl2(I2) + llQ(-> 
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where the constant c depends on u and its derivatives. By reasoning as we did with 
B, we can then get the bound 

|(||p|Il^(^) + IIQIIl^(^))< 

c(l|B(-,0llL2(^) + ||P||^2(^) + ||Q||^2(^)+^^g?+g?+g2j.). 

Via Gronwall's inequaUty and the bounds on ||B||^2 we find 

m-,t)\\l2^n) + m-,t)\\l2(n)<Const, 

where the constant depends on the {Q) norm of Bq and g and on u and its deriva- 
tives. This means that we have an energy estimate 

< (l|Bo||//i(r2) + ll§ll//i(^r2x(o,0)) ' 

where Q is a finite constant depending on u and its derivatives. 

Then, for a general initial data, velocity fields and boundary conditions, we can 
use a standard approximation argument ([ ]) and the above estimate to pass to the 
limit and prove the existence and uniqueness of weak solutions. 

Remark 2.1 The above theorem has been proved in the unit square. It can be easily 
extended to domains with smooth (i.e., boundaries) by using cut-off functions and 
mappings between the domain and the upper-half space. See [ ] and other references 
therein for details. 



3 Semi-discrete Schemes 

As stated before, we will approximate (2.1) with SBP-SAT finite difference schemes. 
We start by defining a SBP operator approximating the first derivative of a continuous 
function w{x) in one space dimension. Let {xi}^^Q be equidistant points in [0, 1] such 
that Xi = ih where h= l/n. We organize the values of w at {xi} in a vector {t) = 
(wo, w„) where Wi = w{xi). Then , we define. 

Definition 3.1 A difference approximation (given by a (^ + 1) x (^ + 1) matrix D) 
for the first derivative is called a Summation- By -Parts (SBP) operator if D = P~^Q for 
^x^matricesPand g, whereP > 0, P = P^ and Q + g^ = ^ = diag(-l,0,0, . . . ,0,0, 1). 

Moreover, P must define a scalar product (w, v) = w^Pv for which the correspond- 
ing norm, ||w||p = (w,w), is equivalent to the standard /^-norm, ||w||2 = /^Lf=i w?. 

SBP operators of different orders of accuracy are presented in several papers, see 
the references in e.g. [ ]. To discretize (2.1), we introduce equidistant meshes in 
the X and y directions with N and M mesh points and Ax = 1/N and Ay = 1/M. 
The discrete solution consists of a column vector of length 2(A/^ + 1) (M + 1) denoted 
y = (y\y^)^, where is a vector of length (A^+ 1)(M+ 1) ordered as 

yi - fyi yi yi yi yi \ 
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and Vf j is the discrete solution at (xi^yj) for ^ = 1, 2. We will use the norm 

where we have introduced the Kronecker product, which is defined as follows: 

Let A and Cbenxn matrices and B and Dbemxm matrices. Then A (8) 5 is the 
nm X nm matrix 

f a\iB . . . a\nB\ 

(A^B) = 

\an\B . . . annBj 

Furthermore, the following rules hold; {A^B){C^D) = (AC^BD), (A (g) 5) + (C (g) 
D) = (A + C) (g) (B^D) and (A 05)^ = (A^ (8)5^). 

To define discrete boundary conditions, we need some further notation. For real 
numbers a/, introduce the 2 x 2 matrices 

where the h is the 2 x 2 identity matrix and the numbers a/ are to be determined later. 
We also need (M+ 1) x (M+ 1) matrices Fq^ and F^j 





/I . . . 








• •1\ 




• • • • 







• • 


• . 




1 • • • 





5 Fnj = 


• • 


• • 1 




• • • • 







• • 


• • 




V 


•J 




I •• 


•7 


and (A^+ 1) X (A^+ 1) matrices F,,o 


and Fjc^M, 








/I • • • 


o\ 




/o.. 


..1\ 




• • • • 







• • 


• • 




1 • • • 





-> Fx,M = 


• • 


. . 1 




• • • • 







• • 


• • 




V 






v- •• 


•7 



Next set 



EjcO = Fro0ly, and E. 



'x,M — ^x,M ^ 



where and ly slvq (N -\- 1) x {N -\- 1) and (M + 1) x (M + 1) identity matrices respec- 
tively and define 



A;, =/2 0diag (4,0^4,1 
Ay=l2 diag {ulo.ul i , . . . , ul^, u\ q, uIm) • 
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Define the matrix C by, 

/- (/. ® (Py-'Qy) ) {h ® {Py-'Qy) ) \ 

\ {{P-'Q,)®Iy)u^ -{{p-'Q,)®Iy)u')- 

Let g be a column vector of the same length as V , where we store the boundary values 
at the appropriate places. Then we can describe the SBP-SAT scheme as 

dtV^A, {h ^{p-^Q^)0ly)V^ Ay {h 0lx^{P~^Qy))V^ CV 

= 1:0,3; ^ (^x"^ ^h) ^^Oy {V - g)^ENy {P-^ 0ly) 0Enj {V - g) (3.1) 

-^E^o0{l^^p-^)0E^o{V-g)^E^^N^{lx^P^^)^E^N{y-g). 



Theorem 3.1 Assume that the velocity field u is a constant given by vi = {u^ ^u^Y , 



and let V be the semi-discrete solution defined by the scheme (3.1). Let u 
and u^~ = {1/ A O),for ^=1,2. If the penalty parameters satisfy 



, (7^ < 

2 ' ^- 2 



and 04 < 



there exists positive constants a and K such that 

t 

||V(Oll'<l|Bo||'+ / / ^{t,x)dxdx, 
J Jdn 



(m^VO) 



(3.2) 



(3.3) 



and the scheme (3.1) is stable. 

Proof The proof is similar to the standard way of proving stability of SBP-SAT 
schemes (see [ ]) and follows the proof for obtaining energy stability of the con- 
tinuous problem in theorem 2.1. We outline the proof for the sake of completeness. 
For simplicity, we consider the case of constant velocities by setting C = in (3.1) 
We start by multiplying (3.1) with V^{l20Px^Py) to obtain. 



{l2^Px^Py)dtV 



V {A,0Q,0Py)V^V' {Ay( 



^Qy)y 



-V^ih^Px^Py)^ 



{Eoy ^P^^^Iy) Eoj + {En J (8) ^ ( 



)Iy) En J 



^{E^^O^I^^Py ')E^^o^{E^^N^Ix^Py 'R,M 

Adding this to its transpose and using the definition of SEP operators, we obtain 

I|2 



V. 
(3.4) 



dt 



V 



-V^ (Ai< 
-2V^{l2 



)Py)V^V^ {A2^Px^^y)V 



)Py)i 



{EOJ^P-^ (g)/y) Eoj + {Enj^P~^ ^ly) EN^y 
+ {E^^o ^I^^P-^) E^^o + {E^^N 0lx^Py~^) ^x,A 



V, 
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which impHes 




+ {V^fl) ^ Px (V^fi) - (vIn) ^ Px {V^,n) 

+ 2 [^1 ( (^0.^) ^ (^Oy) + {yO,y) ^ Py Ky) ) + ^2 {V^y) ^ Py (4,) 



Using (3.2) and integrating in time gives the energy estimate (3.3). 

Remark 3.1 The above proof of stabiHty assumes a constant velocity field. A proof 
of stability with a general velocity fields has been obtained in a recent paper [ ' ' ] by 
using the principle of frozen coefficients. The resulting stability estimate will lead to 
an exponential growth of energy (similar to (2.6)) due to the presence of lower order 
terms. 

We conclude this section with a few comments. For simplicity, we have only con- 
sidered Cartesian meshes. However, the proofs are readily generalized to curvilinear 
grids by transforming the domain to a Cartesian. A stability proof is then obtained by 
freezing the coefficients. However, that requires P to be diagonal, [ ]. Furthermore, 
multi-block grids can also be handled and stable interfaces derived in a similar way 
as in, [l o]. 

4 Numerical Experiments 

We test the SBP-SAT schemes of the previous section on a suite of numerical exper- 
iments in order to demonstrate the effectiveness of these schemes. We will use two 
different schemes : SBP2 and SBPA scheme which are second-order (first-order) and 
fourth order (second-order) accurate in the interior (boundary) resulting in an overall 
second and third-order of accuracy. Time integration is performed by using a second 
order accurate Runge-Kutta scheme at a CFL number of 0.45 for all numerical ex- 
periments. We found that using a fourth order accurate Runge-Kutta scheme resulted 
in negligible differences in the numerical results. The schemes have bounded errors, 
a typical behavior for hyperbolic equations with characteristic boundary conditions 
as shown in [ ] . Errors are propagated through the domain and leave the domain on 
account of the transparent boundaries. Hence, errors do not accumulate in time. On 
small domains, spatial errors become dominant. 

Numerical experiment 1: In this experiment, we consider (2.1) with the divergence- 
free velocity field vi{x^y) = (— The exact solution can be easily calculated by 
the method of characteristics and takes the form 



B{x,t)=R{t)Bo{R{-t)x), 



(4.1) 
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where R{t) is a. rotation matrix with angle t and represents rotation of the initial data 
about the origin. 

We consider the same test setup as in [^^] and [v.] by choosing the divergence free 
initial data, 

Bo(x,y)=4(^7,).-2o((-i/2)^+/), (4.2) 

and the computational domain [—1,1] x [—1,1]. Since the exact solution is known in 
this case, one can in principle use this to specify the boundary data g. Instead, we 
decided to mimic a free space boundary (artificial boundary) by taking g = 0. (which 
is a good guess at a far-field boundary). 

We run this test case with SBP2 and SBP4 schemes and present different sets of 
results. In Figure 4.1, we plot |B| = {\B^ p + ^2 at times t = 71 (half-rotation) 
and t = 271 (one full rotation) with the SBP2 and SBP4 schemes. As shown in this 





(a) half rotation, SBP2 (b) full rotation, SBP2 





I: 

i: 

(c) half rotation, SBP4 (d) full rotation, SBP4 

Fig. 4.1 Numerical results for |B| in experiment 1. 



figure, SBP2 and SBP4 schemes resolve the solution quite well. In fact, SBP4 is very 
accurate and keeps the hump intact throughout the rotation. In Table 4.1, we present 
percentage relative errors in The errors are computed at time t = 271 (one rotation) 
on a sequence of meshes for both the SBP2 and SBP4 schemes. The results show that 
the errors are quite low, particularly for SBP4 and the rate of convergence approaches 
the expected values of 2 for SBP2 and 3 for SBP4. Furthermore, the order of accuracy 
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Grid size 


SBP2 


rate 


SBP4 


rate 


40x40 


6.9-101 




8.0 -10^ 




80x80 


2.1-101 


1.7 


5.0-10-1 


4.0 


160x160 


5.5 -10^ 


2.0 


4.5-10-2 


3.5 


320x320 


1.3 -10^ 


2.0 


5.1-10-3 


3.1 


640x640 


3.3-10-1 


2.0 


6.4-10-4 


3.0 



Table 4.1 Relative percentage errors in for |B| at time t ^2n and rates of convergence for numerical 
experiment 1 with SBP2 and SEP A schemes. 

is unaffected at these resolutions by using zero Dirichlet boundary data instead of the 
exact solution at the boundary. 

In order to compare the SBP schemes of this paper with other existing schemes, 
we choose to compute the solutions for this problem with both the first- and the 
second-order divergence-preserving scheme of [ ], which we label as the TF and 
TF2 schemes. Furthermore, we compute the solutions using the first order stable 
upwind scheme designed in [ ], labeled the SUS scheme. The relative errors with 
each of these schemes are shown in Table 4.2. Results in Tables 4.1 and 4.2 show 



Grid size 


SUS 


TF 


TF2 


40x40 


8.6-101 


7.6-101 


1.8-101 


80x80 


7.3-101 


6.4-101 


1.3-101 


160x160 


5.4-101 


4.7-101 


3.0- 10^ 


320x320 


3.6-101 


3.3-101 


1.0-10^ 


640x640 


2.0-101 


1.4-101 


2.7-10-1 



Table 4.2 Relative percentage errors in for |B| at f = 2ti and for numerical experiment 1 with the SUS, 
TF, TF2, SBP2 and the SBPA schemes. 



that the TF and SU S schemes lead to similar errors and these errors are considerably 
larger than the errors generated by the TF2 and SBP2 schemes, while the errors 
generated by the SBP4 scheme are much smaller again. 

A fair comparison of the the five schemes SUS, TF, TF2, SBP2 and SBP4 requires 
information on the computational work with each scheme for the same error level. 
We observe from tables 4.1 and 4.2 that for a given relative error of approximately 
20 percent, the first-order SUS scheme requires a 640 x 640 mesh, the TF scheme 
requires a 500 x 500 mesh (based on extrapolation from table 4.2), whereas both 
the second-order schemes require meshes coarser than a 50 x 50 mesh. The fourth- 
order scheme yields similar error levels on even coarser meshes. Thus, the second- 
order schemes require about 1 % of the total grid points to the first-order schemes to 
produce comparable errors. Even taking into account that the second order schemes 
use more operations for each grid point, this still makes the second order schemes 
at least 25 — 30 times more efficient than the first order schemes. Similarly an error 
level of about one percent is attained with SBP2 on a 320 x 320 mesh, with TF2 
on a similar 320 x 320 mesh and with SBP4 on a 50 x 50 mesh. Thus the second 
order schemes need about 36 times more grid points to produce errors similar to 
those of the fourth order schemes. Taking extra work for the fourth-order scheme per 
grid point into account, we still get that the fourth-order scheme is roughly 10 times 
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more efficient than the second-order schemes. These numbers are approximations 
but display a clear qualitative trend i.e., it is much more efficient to use high-order 
schemes for the induction equations. 

As the solution (4. 1) in this case is smooth, it is also a solution for the constrained 
form (1.3). Furthermore, the initial data is divergence free and so is the exact solu- 
tion. We did not attempt to preserve any particular form of discrete divergence while 
designing the SBP schemes (3.1). A natural thing would be show that some form of 
discrete divergence produced by the schemes was bounded in We were unable to 
obtain such a divergence bound for (3. 1) in this paper. A related SBP-SAT scheme for 
the "conservative" symmetric form (1.5) with SBP operators for discretizing spatial 
derivatives coupled with a novel discretization of the source terms in (1 .5) was shown 
to have bounded discrete divergence in a recent paper [ ]. 

In the absence of a rigorous divergence bound, we proceed to examine how diver- 
gence errors generated by the SBP schemes behave and whether they have any impact 
on the quality of the discretization. We define the following discrete divergence, 

divp(y) = {p-'Q,®iy)y' ^{h®p-'Qy)v\ 

This corresponds to the standard centered discrete divergence operator at the corre- 
sponding order of accuracy. The divergence errors in and rates of convergence at 
time t = 271 for the SBP2 and SBP4 schemes on a sequence of meshes are presented 
in Table 4.3. From Table 4.3, we conclude that although the initial divergence is zero. 



Grid size 


SBP2 


rate 


SBP4 


rate 


20x20 


1.0-10" 




7.3-10-1 




40x40 


8.0-10-1 


0.4 


1.2-10-1 


2.6 


80x80 


2.7-10-1 


1.6 


8.2-10-3 


3.8 


160x160 


7.0-10-2 


2.0 


1.0-10-3 


3.0 


320x320 


2.5-10-2 


1.5 


1.7-10-4 


2.6 



Table 4.3 Numerical Experiment 1 : Divergence (errors) in and rates of convergence at time t ^2n for 
both the SBP2 and SBPA schemes. 



the discrete divergence computed with both the SBP2 and SBP A schemes is not zero. 
However, the divergence errors are very small even on fairly coarse meshes and con- 
verge to zero at a rate of 1.5 and 2.5 for SBP2 and SBP A schemes respectively. A sim- 
ple truncation error analysis suggests that these rates for the SBP2 and SBPA schemes 
are optimal. The quality of the approximations is good and the rates of convergence 
do not seem to suffer from not preserving any form of discrete divergence. 

In order to compare with existing schemes, we compare the divergence errors gen- 
erated by the SU S, TF and the TF2 schemes with the SBP2 and the SBP4 schemes in 
table 4.4. From Table 4.4, we can draw the following conclusions about divergence 
errors. The SUS scheme is not tailored to preserve any form of discrete divergence. 
The divergence errors generated by this scheme seems to be low on coarse meshes. 
The TF and TF2 schemes are designed to preserve a special form of discrete diver- 
gence which is different from the standard central form. Nevertheless, the analysis 
presented in [ ] suggested that the errors in the standard divergence operator will 
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Grid size 


SUS 


TF 


TF2 




40x40 


l.MO-i 


2.7-10-2 


1.2-10" 


-2 


80x80 


1.3-10-1 


1.7-10-2 


4.0-10- 


-3 


160x160 


1.4-10-1 


1.4-10-2 


2.4-10" 


-3 


320x320 


1.1-10-1 


1.2-10-2 


9.7-10- 


-4 



Table 4.4 Numerical Experiment 1 : The discrete divergence divp in /2 at f = 27i for the SUS, TF and TF2 
schemes. 



also be quite low. This is indeed the case. On the coarser meshes, the divergence is 
much larger for the SBP2 scheme than the TF schemes, but from Table 4.2 we see 
that the errors in the solution are similar. 

Furthermore, the divergence errors converge quickly for the SBP4 scheme, as 
well as as the for the TF2 scheme. The above results indicate that controlling some 
form of discrete divergence is not necessary to approximate solutions of the magnetic 
induction equations in a stable and accurate manner. 

Next, we consider long time integration. The energy estimate (3.3) suggests that 
the energy of the approximate solutions can grow exponentially in time. In order 
to test this we computed approximate solutions with the SBP2, SBP4 and the TF2 
schemes till time t = 100;r, i.e., for fifty full rotations on a 100 x 100 mesh. The 
numerical results in are presented in Figure 4.2 and Table 4.5. These computations 



© 











(a) r = (b) t = IOtt, SBP2 (c) t = IOOtt, SBP4 

Fig. 4.2 Numerical results for |B| in experiment 1. 



27lt 


SBP2 


SBPAr 


TF2 


t^l 


2.1-101 


5.1-10-1 


8.8-10" 




7.7-101 


2.7 - 10^ 


3.2-101 


t^lO 


1.0-102 


4.7 - 10^ 


5.0-101 


t^l5 


1.1-102 


6.6 -10^ 


6.3-101 


t^20 


1.2-102 


8.7 - 10^ 


7.2-101 


t^30 


1.2-102 


1.9-101 


8.4-101 


t^40 


1.3-102 


3.1-101 


9.2-101 


f = 50 


1.4-102 


4.3-101 


1.0-102 



Table 4.5 Relative percentage /2 errors in |B| with SBP2, SEP A and TF2 for numerical experiment 1. 



were performed on a fixed 100 x 100 mesh. In Figure 4.2, we compare the SBP2 and 
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SBP4 schemes after five and fifty rotations respectively. We see that after 5 rotations, 
SBP2 gives a "hump" which is somewhat smeared and with a pronounced asymmetry. 
On the other hand, the hump produced by the SBP4 scheme is much more accurate. 
As shown in Table 4.5, the absolute errors with the SBP4 scheme are much lower 
than the errors due to the second-order schemes SBP2 and TF2. In fact, the errors 
with SBP2 after just five rotations are about three times the error with SBP4 after 
fifty rotations. This experiment makes a strong case for using high-order schemes for 
problems requiring long time integration. 

Numerical Experiment 2: In the previous numerical experiment, the hump was con- 
fined to the interior of the domain during the rotation. Hence, the choice of zero 
Dirichlet data at the boundary was reasonable and led to stable and accurate approxi- 
mations. In order to illustrate the effect of the boundary better, we choose the compu- 
tational domain [0, 1] x [0, 1] and use the same velocity field and initial data as in the 
previous experiment. Now, the hump "exits" the domain at one part of the boundary 
(including a corner) and will re-enter the domain from another part of the boundary. 
The choice of boundary discretization becomes crucial in this case. 

We select the exact solution (4.1) restricted to the boundary as the Dirichlet 
boundary data in (3.1). In Figure 4.3, the approximate solutions computed with both 





(c) SBP4,r = 71 /2 



Fig. 4.3 Numerical results for experiment 2. Mesh size 100 x 100. 



(d) SBP4,r = 2k 



SBP2 and SBP4 on a 100 x 100 mesh are plotted at time t = n/2 (quarter rotation) 
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and time t = 271 (full rotation). As shown in this figure, both schemes perform very 
well. The hump at both the exit as well as the re-entry is clearly resolved with no 
noticeable numerical artefacts or reflections. 



Grid size 


SBP2 


rate 


SBP4 


rate 


10x10 


2.5 -10^ 




1.1 -10^ 




20x20 


5.8 -10^ 


2.1 


1.5 -10° 


2.9 


40x40 


1.3-10^ 


2.0 


1.6-10-1 


3.3 


80x80 


3.0-10-1 


2.0 


1.6-10-2 


3.2 


160x160 


7.4-10-2 


2.0 


1.9-10-3 


3.1 



Table 4.6 Numerical experiment 2: Relative percentage errors for |B| in /2 and rates of convergence for 
both SBP2 and SBP4. 



Grid size 


SBP2 


rate 


SBP4 




rate 


10x10 


6.4-10-1 




9.7-10" 


-2 




20x20 


3.9-10-1 


0.7 


2.4-10" 


-2 


2.0 


40x40 


9.1-10-2 


2.2 


1.9-10" 


-3 


3.6 


80x80 


2.6-10-2 


1.8 


3.0-10" 


-4 


2.7 


160x160 


8.9-10-3 


1.6 


5.1-10" 


-5 


2.5 



Table 4.7 Numerical experiment 2: Divergence (errors) in /2 and rates of convergence for both SBP2 and 
SBP4 at time t = 2k. 



As shown in Table 4.6, the errors are low after one full rotation for both the SBP2 
and SBP4 schemes. In fact, the size of relative errors is lower than in the previous 
numerical experiment. As expected, the rates of convergence tend to 2 and 3 for SBP2 
and SBP4 respectively. In Table 4.7 the divergence errors and their convergence rates 
are listed. They are small and the convergences approach the expected values 1.5 and 
2.5. 

On the other hand, when we tried to compute this example with the divergence 
preserving TF and TF2 schemes, the solution blew up on account of boundary insta- 
bilities. 

Numerical Experiment 3: (Discontinuous solutions.) As remarked earlier, the mag- 
netic induction equations (2.1) are a sub-model in the nonlinear MHD equations. As 
a consequence, one must solve the induction equation with both discontinuous ve- 
locity fields and initial data. It is therefore interesting to see how well the SBP-SAT 
schemes handle discontinuous velocity fields and initial data. 

The SBP operators use centered finite differences in the interior. It is well known 
that using central differences leads to oscillations around discontinuities. Therefore 
the SBP schemes cannot be used directly in this regime, see [12] for details. To cal- 
culate solutions with discontinuities, one adds a small amount of explicit numerical 
diffusion that retain the accuracy of the first derivative SBP approximations as well as 
maintain the energy stability of the SBP scheme. We will use these operators together 
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with the SBP2 and SBPA schemes in order to compute discontinuous solutions of the 
magnetic induction equations. 

The second-order (fourth-order) SBP operator for the first derivative with a second- 
order (fourth-order) numerical diffusion operator gives an approximation which is 
formally second-order (fourth-order) accurate in the interior of the computational 
domain. It turns out that a different scaling (dividing by the mesh size) of the numer- 
ical diffusion operator leads to a first order (third-order) "upwind" scheme. We will 
test all these numerical diffusion operators a numerical experiment first described in 



The computational domain is [0, 1] x [0, 1]. Consider the constant velocity field, 
u = (1,2)^ and the discontinuous initial data, 



The initial discontinuity simply moves along the diagonal of the domain. We use the 
exact solution restricted to the boundary as the Dirichlet boundary data. Tests with 
generic SBP -SAT schemes, (3.1), showed that the approximate solutions were very 
oscillatory, and we damp these oscillations by adding numerical diffusion. 

We test the SBP2 (SBP4) scheme with the standard second-order (fourth-order) 
numerical diffusion operator as well as the scaled numerical diffusion operator to 
obtain the first-order (third-order) SB PI and SBP3 schemes. The results on a 100 x 
100 mesh at time t = 0.5 are plotted in Figure 4.4. A plot at this time is of interest 
as some part of the solution has interacted with the boundary and exited the domain, 
whereas most of the front is still inside the domain. From Figure 4.4, we see that the 
boundary discretization works well in all cases and does not lead to any significant 
oscillations in the domain. The SB PI scheme is the most dissipative with significant 
smearing at the discontinuity. However, this scheme also has no over/under shoots or 
oscillations and the solution is TVD. The SBP2 scheme with second-order numerical 
diffusion operator is oscillatory near the discontinuity with dispersive waves on both 
sides of it. The smearing is considerably less than that of the SBPl scheme. The SBP4 
scheme with standard fourth-order numerical diffusion is even more oscillatory and 
leads to a larger overshoot. The SBP3 scheme damps these oscillations somewhat and 
still keeps the sharpness at the discontinuity making it an acceptable alternative. 



5 Conclusion 

We have considered the magnetic induction equations that arise as a submodel in 
the MHD equations of plasma physics. Various forms of these equations were pre- 
sented including the symmetric forms that are well-posed with general initial data 
and Dirichlet boundary conditions. 



[6]. 




In this case, the exact solution (see [ ]) of (2.1) reads 



B{x,y,t)=Bo{x-t,y-2t). 
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(a) SBPl (b) SBP2, second-order diffusion 




(c) SBP3 (d) SBP4, fourth-order diffusion 

Fig. 4.4 Numerical results for B^ {x,y, 0.5) in experiment 3. 

Standard numerical methods of the finite difference/finite volume type have dealt 
with discretizations of the constrained form (1.3) and attempted to preserve some 
form of the divergence constraint. 

We describe SBP-SAT based finite difference schemes for the initial- boundary 
value problem corresponding to the magnetic induction equations. These schemes 
were based on the non-conservative symmetric form (1.4) and use SBP finite differ- 
ence operators to approximate spatial derivatives and a SAT technique for implement- 
ing boundary conditions. The resulting schemes were energy stable and higher order 
accurate. 

These schemes were tested on a series of numerical experiments, which illus- 
trated their stability and high-order of accuracy. Interesting solution features were 
resolved very well. The fourth-order scheme was found to be well suited for long 
time integration problems. Despite the fact that the schemes were not preserving any 
particular form of discrete divergence as well as the lack of a rigorous discrete di- 
vergence bound, the divergence errors generated by the schemes were quite low and 
converged to zero at the expected rates when the mesh was refined. The schemes 
were compared with two existing lower order schemes and one divergence preserv- 
ing second order scheme. Despite lacking any divergence bounds, the SBP schemes 
performed at least as well as the schemes with a divergence bound. 

The numerical experiments indicate that the SBP-SAT framework is effective in 
approximating solutions of the magnetic induction equations to a high order of accu- 
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racy. In the future we plan to extend these schemes to magnetic induction equations 
with resistivity. 

References 

1. D.S. Balsara and D. Spicer. A staggered mesh algorithm using high order Godunov fluxes to ensure 
solenoidal magnetic fields in magnetohydrodynamic simulations. / Comp. Phys., 149(2): 270-292, 
1999. 

2. N. Besse and D. Kroner. Convergence of the locally divergence free discontinuous Galerkin methods 
for induction equations for the 2D-MHD system. M2 AN Math. Model. Num. Anal 39(6): 1177-1202, 
2005. 

3. J.U. Brackbill and D.C. Barnes. The effect of nonzero div5 on the numerical solution of the magne- 
tohydrodynamic equations. /. Comp. Phys., 35:426-430, 1980. 

4. W. Dai and P.R. Woodward. A simple finite difference scheme for multi-dimensional magnetohydro- 
dynamic equations. /. Comp. Phys., 142(2):331-369, 1998. 

5. C. Evans and J.F. Hawley. Simulation of magnetohydrodynamic flow: a constrained transport method. 
Astrophys. J., 332:659, 1998. 

6. F. Fuchs, K.H. Karlsen, S. Mishra and N.H. Risebro. Stable upwind schemes for the Magnetic Induc- 
tion equation. Preprint, Submitted. 

7. F. Fuchs, S. Mishra and N.H. Risebro. Spfitting based finite volume schemes for ideal MHD equations. 
in press Jl. Comput. Phys., 

8. S.K. Godunov. The symmetric form of magnetohydrodynamics equation. Num. Meth. Mech. Cont. 
Media, 1:26-34, 1972. 

9. B. Gustafsson, H.-O. Kreiss, and J. Oliger. Time dependent problems and difference methods. John 
Wiley & Sons, Inc., 1995. 

10. H.O. Kreiss and J. Lorenz. Initial-Boundary value problems and the Navier-Stokes equations. Aca- 
demic Press, Boston, 1989. 

11. R.J. LeVeque. Finite volume methods for hyperbolic problems. Cambridge university press, Cam- 
bridge, 2002. 

12. K. Mattsson, M. Svard and J. Nordstrom. Stable and Accurate Artificial Dissipation. Journal of 
Scientific Computing, Vol.21,No.l, August 2004. 

13. K. Mattsson and J. Nordstrom. Summation by parts operators for finite difference approximations of 
second derivative. Journal of Computational Physics, 199(2004), 503-540. 

14. S. Mishra and M. Svard. On stability of numerical schemes via frozen coefficients and magnetic 
induction equations. Preprint, Submitted. 

15. J. Nordstrom. Error Bounded Schemes for Time-dependent Hyperbolic Problems SI AM J. Sci. Corn- 
put., 30(2^1), 46-59. 

16. J. Nordstrom and M. H. Carpenter High-Order Finite Difference Methods, Multidimensional Linear 
Problems, and Curvilinear Coordinates /. Comput. Phys, 173(2001), 149-174. 

17. K.G. Powell. An approximate Riemann solver for magneto-hydro dynamics (that works in more than 
one space dimension). Technical report, 94 -24, ICASE, Langley, VA, 1994. 

18. K.G. Powell, P.L. Roe. T.J. Linde, T.I. Gombosi and D.L. De Zeeuw, A solution adaptive upwind 
scheme for ideal MHD. J. Comp. Phys, 154(2), 284 - 309, 1999 

19. J. Rauch. Partial Differential Equations, Springer, 1991. 

20. D.S. Ryu, F. Miniati, T.W Jones and A. Frank. A divergence free upwind code for multidimensional 
magnetohydrodynamic flows. Astrophys. J., 509(l):244-255, 1998. 

21. M. Svard On coordinate transformations for summation-by-parts operators /. Sci. Comput. 20(2004), 
29-42. 

22. M. Svard and J. Nordstrom. On the order of accuracy for difference approximations of initial- 
boundary value problems. Journal of Computational Physics, 218(2006), 333-352. 

23. M. Torrilhon and M. Fey. Constraint-preserving upwind methods for multidimensional advection 
equations. SIAM. J Num. Anal, 42(4): 1694-1728, 2004. 

24. G. Toth. The div5 = constraint in shock capturing magnetohydrodynamics codes. /. Comp. 
P/iy^., 161:605-652, 2000. 



